      SUBROUTINE addmult(a,b,c,d,n)
      integer n
      double precision a(n), b(n), c(n), d(n)
      integer i
      DO i=1, n
        a(i) = b(i) + 1.3*c(i) + 1.7*d(i)
      END DO
      RETURN
      END

      SUBROUTINE do_rk4(a,b,c,d,d1,d2,n)
      integer n
      double precision a(n), b(n), c(n), d(n)
      double precision d1,d2
      integer i
      DO i=1, n
        a(i)=a(i)+d1*d(i)
        b(i)=c(i)+d2*d(i)
      END DO
      RETURN
      END

      !always c storage
      SUBROUTINE vxyz3(a,b,n1ind,n1,n2ind,n2,n3ind,n3,ntot)
      integer n1,n2,n3,ntot
      integer n1ind(n1),n2ind(n2),n3ind(n3)
      complex*16 a(ntot), b(ntot)
      complex*16 ig,jg,kg

      !internal
      integer i,j,k,ijk

      ijk=1
      DO i=1,n1
        ig=cmplx(0,1.1*n1ind(i))
        DO j=1,n2
          jg=cmplx(0,1.1*n2ind(j))
          DO k=1,n3
            kg=cmplx(0,1.1*n3ind(k))
            a(ijk  )=jg*b(ijk+2)-kg*b(ijk+1)
            a(ijk+1)=kg*b(ijk  )-ig*b(ijk+2)
            a(ijk+2)=ig*b(ijk+1)-jg*b(ijk  )
            ijk=ijk+3
          ENDDO
        ENDDO
      END DO
      RETURN
      END

      !assume c storage
      SUBROUTINE vxy3z(a,b,n1ind,n1,n2ind,n2,n3ind,n3)
      integer n1,n2,n3
      integer n1ind(n1),n2ind(n2),n3ind(n3)
      complex*16 a(n3,3,n2,n1), b(n3,3,n2,n1)
      complex*16 ig,jg,kg
      integer i,j,k

      DO i=1,n1
        ig=cmplx(0,1.1*n1ind(i))
        DO j=1,n2
          jg=cmplx(0,1.1*n2ind(j))
          DO k=1,n3
            kg=cmplx(0,1.1*n3ind(k))
            a(k,1,j,i)=jg*b(k,3,j,i)-kg*b(k,2,j,i)
            a(k,2,j,i)=kg*b(k,1,j,i)-ig*b(k,3,j,i)
            a(k,3,j,i)=ig*b(k,2,j,i)-jg*b(k,1,j,i)
          ENDDO
          !DO k=1,n3
          !  kg=cmplx(0,1.1*(k-1))
          !  a(k,2,j,i)=kg*b(k,1,j,i)-ig*b(k,3,j,i)
          !ENDDO
          !DO k=1,n3
          !  a(k,3,j,i)=ig*b(k,2,j,i)-jg*b(k,1,j,i)
          !ENDDO
        ENDDO
      END DO
      RETURN
      END

      SUBROUTINE v3xyz(a,b,n1ind,n1,n2ind,n2,n3ind,n3)
      integer n1,n2,n3
      integer n1ind(n1),n2ind(n2),n3ind(n3)
      complex*16 a(n3,n2,n1,3), b(n3,n2,n1,3)
      complex*16 ig,jg,kg
      integer i,j,k,d

      DO i=1,n1
        ig=cmplx(0,1.1*n1ind(i))
        DO j=1,n2
          jg=cmplx(0,1.1*n2ind(j))
          DO k=1,n3
            kg=cmplx(0,1.1*n3ind(k))
            a(k,j,i,1)=jg*b(k,j,i,3)-kg*b(k,j,i,2)
          ENDDO
        ENDDO
      END DO
      DO i=1,n1
        ig=cmplx(0,1.1*n1ind(i))
        DO j=1,n2
          jg=cmplx(0,1.1*n2ind(j))
          DO k=1,n3
            kg=cmplx(0,1.1*n3ind(k))
            a(k,j,i,2)=kg*b(k,j,i,1)-ig*b(k,j,i,3)
          ENDDO
        ENDDO
      END DO
      DO i=1,n1
        ig=cmplx(0,1.1*n1ind(i))
        DO j=1,n2
          jg=cmplx(0,1.1*n2ind(j))
          DO k=1,n3
            kg=cmplx(0,1.1*n3ind(k))
            a(k,j,i,3)=ig*b(k,j,i,2)-jg*b(k,j,i,1)
          ENDDO
        ENDDO
      END DO
      RETURN
      END
